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Abstract 

This  paper  addresses  the  development  of  a  temperature-dependent  constitutive  model  for  relaxor  fer¬ 
roelectric  materials.  These  compounds  exhibit  a  diffuse  transition  region  between  the  paraelectric  and 
ferroelectric  phases  due  to  the  chemical  heterogeneity  of  the  materials.  At  low  temperatures,  the  materials 
exhibit  significant  dielectric  hysteresis  in  the  relation  between  the  applied  field  E  and  the  macroscopic 
polarization  P,  with  the  degree  of  hysteresis  decreasing  as  the  temperature  increases  to  the  freezing  tem¬ 
perature  Tf.  Above  the  freezing  temperature,  the  relation  is  single- valued  but  nonlinear. 

These  phenomena  are  modeled  by  assuming  that  the  material  is  comprised  of  an  aggregate  of  microp- 
olar  regions  having  a  range  of  Curie  temperatures.  Thermodynamic  principles  are  employed  to  obtain  a 
micropolar  model  which  predicts  the  saturation  polarization  and  distribution  of  regions  as  a  function  of 
temperature.  A  corresponding  macroscopic  model  is  then  constructed  to  predict  the  dielectric  behavior  of 
the  material  above  the  freezing  temperature.  Hysteresis  below  the  freezing  point  is  incorporated  through 
the  quantification  of  energy  required  to  bend  and  translate  domain  walls  pinned  at  inclusions  inherent  to 
the  material.  The  resulting  ODE  model  quantifies  the  constitutive  nonlinearities  and  hysteresis  exhibited 
by  the  materials  through  a  wide  range  of  temperatures  and  input  drive  levels.  The  predictive  capabil¬ 
ities  of  the  model  are  illustrated  through  a  comparison  with  PMN  and  PMN-PT-BT  data  collected  at 
temperatures  ranging  from  133°  K  to  313°  K. 
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1  Introduction 


Relaxor  ferroelectric  compounds  are  proving  advantageous  in  a  number  of  high  performance  applications 
due  to  their  capability  for  generating  large  strains  with  reduced  hysteresis  when  employed  at  temperatures 
within  their  diffuse  transition  region.  For  example,  a  lead  magnesium  composition  Pb(Mg1//3,  Nli^/.-; ) O3 
(PMN-PT)  has  been  demonstrated  to  produce  at  least  a  6  dB  increase  in  acoustic  source  levels  over  Pb(Zr, 
Ti)03  (PZT)  when  employed  as  drivers  for  underwater  sonar  transducers  [12,  15].  The  capability  of  PMN 
for  applications  which  require  micropositioning  is  illustrated  by  its  success  in  adaptive  optic  systems.  When 
employed  in  multi-layered  arrays  for  deformable  mirrors,  PMN  has  increased  mirror  sensitivity  from  0.3 
nrn/V  to  12.0  nm/V  relative  to  monolithic  PZT  [6,  9,  10]. 

The  stoichiometries  and  electromechanical  mechanisms  which  provide  the  materials  with  superior  per¬ 
formance  capabilities  also  imbue  them  with  a  complex  constitutive  behavior  which  must  be  quantified  to 
achieve  their  potential  in  high  performance  system  design.  Throughout  the  temperature  range,  the  mate¬ 
rials  exhibit  nonlinear  saturation  effects  and  positive  induced  strains  due  to  their  electrostrictive  nature. 
Hence  they  must  be  biased  through  an  applied  field  so  that  they  operate  in  one  quadrant  of  the  field-strain 
curve.  The  incorporation  of  the  inherent  nonlinearities  provides  the  capability  for  designing  tunable  sensors 
once  the  nonlinearity  is  characterized.  Additionally,  the  materials  exhibit  temperature-dependent  hystere¬ 
sis  below  the  freezing  temperature  with  the  degree  of  hysteresis  increasing  as  temperatures  decrease  as 
illustrated  in  Figure  1.  Finally,  the  strains  generated  by  the  materials  are  dependent  on  both  the  operating 
temperature  and  frequency. 

The  manner  through  which  these  properties  are  manifested  in  transducer  design  can  be  illustrated 
through  the  consideration  of  a  sonar  projector.  When  driven  at  high  fields  and  high  frequency,  PMN  will 
heat  due  to  dielectric  hysteresis.  Since  PMN’s  hysteresis  decreases  with  temperature,  actuators  reach  an 
equilibrium  temperature  when  heat  dissipation  due  to  conduction,  convection  and  radiation  balances  the 
internal  heat  generated  by  the  ceramic.  The  final  equilibrium  temperature  depends  on  the  configuration 
of  the  actuator,  its  electric  drive  cycle,  and  its  surrounding  environment.  Shankar  and  Horn  [16]  predicted 
a  40°  C  rise  for  a  flextensional  sonar  transducer  submersed  in  water  using  a  simple  heat  generation  model 
for  PMN.  These  predictions  agreed  quantitatively  with  experimental  measurements  made  on  an  actual 
transducer. 

To  minimize  the  deleterious  effects  of  self  heating  and  to  optimize  the  design  capabilities  of  transducers 
which  utilize  relaxor  ferroelectrics,  it  is  necessary  to  develop  temperature-dependent  constitutive  relations 
which  quantify  the  nonlinear  behavior  and  hysteresis  throughout  the  operating  range  of  the  materials. 
The  initial  models  proposed  by  Horn  and  Shankar  [11]  quantified  the  electrostrictive  strains  and  nonlinear 
saturation  behavior  of  the  materials  for  isothermal  conditions  and  transition  temperatures  where  hysteresis 
is  negligible.  Subsequent  experiments  by  Brown  et  al.  [2]  demonstrated  that  the  model  accurately  pre- 


Figure  1.  Hysteresis  and  saturation  polarization  of  PMN  as  a  function  of  increasing  temperature  (from 
Glazounov  et  al.  [7]).  Abscissas:  electric  field  (MV/m),  ordinates:  polarization  (C/m2). 
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diets  the  coupled  electromechanical  behavior  of  the  materials  for  isothermal,  hysteresis-free  (anhysteretic) 
operating  conditions.  The  hysteresis  inherent  to  the  materials  at  lower  temperatures,  where  the  material 
exhibits  a  well-defined  domain  structure,  was  then  modeled  by  Smith  and  Horn  [19,  20]  through  the  quan¬ 
tification  of  losses  due  to  domain  wall  motion.  This  provided  a  unified  framework  which  quantified  both 
the  nonlinear  behavior  and  hysteresis  exhibited  by  the  materials  at  fixed  temperatures  and  quasistatic 
operating  conditions.  An  initial  model  characterizing  the  temperature-dependence  of  PMN  was  developed 
by  Glazounov  at  al.  [7]  using  as  a  basis  the  same  polarization  model  employed  by  Horn  and  Shankar. 
This  model  focused  on  the  anhysteretic  behavior  of  the  materials  and  was  based  on  the  assumption  that 
the  material  consists  of  micropolar  regions  having  a  distribution  of  local  Curie  temperatures.  The  model 
quantified  the  saturation  polarization  of  the  aggregate,  the  density  of  micropolar  regions  and  the  distribu¬ 
tion  of  transition  temperatures  for  PMN  but  did  not  account  for  interactions  between  polar  regions.  The 
theory  of  Glazounov  et  al.  was  subsequently  extended  by  Horn  and  Shankar  [13]  to  provide  a  quasistatic, 
temperature-dependent  model  for  PMN  which  quantified  the  nonlinear  behavior  of  the  material  but  did 
not  accommodate  the  hysteresis  present  at  temperatures  below  the  freezing  temperature. 

In  this  paper,  we  extend  the  previous  theory  to  provide  a  unified  model  which  quantifies  the  tem¬ 
perature-dependent  constitutive  nonlinearities  and  hysteresis  exhibited  by  relaxor  ferroelectrics  from  low 
temperature,  ferroelectric  regimes  through  the  diffuse  transition  region  preceding  their  paraelectric  state. 
Like  Glazounov  et  al.,  we  assume  that  the  relaxor  material  consists  of  superparaelectric,  micropolar  re¬ 
gions  with  a  diffuse  spectrum  of  Curie  temperatures.  However,  we  use  a  simple  Ising  model  with  near 
neighbor  ion  interaction  to  represent  the  thermodynamics  of  the  individual  micropolar  regions.  Based  on 
the  “random-layer”  model  of  spatial  B-site  ordering  for  relaxor  ferroelectrics,  we  calculate  the  dispersion 
of  local  Curie  temperatures.  The  micropolar  model  predicts  the  saturation  polarization  and  the  density  of 
the  polar  regions  as  a  function  of  temperature.  A  macroscopic  model,  which  includes  interaction  between 
the  micropolar  regions,  is  then  constructed  to  predict  the  anhysteretic  polarization  behavior  of  the  mate¬ 
rials  above  the  freezing  temperature.  The  hysteresis  inherent  to  the  materials  below  this  temperature  is 
incorporated  through  the  extension  of  the  domain  theory  in  [20]  to  incorporate  thermal  dependence  in  the 
energy  required  to  translate  domain  walls.  As  demonstrated  by  comparisons  to  experimental  PMN  and 
PMN-PT-BT  data,  the  resulting  model  accurately  quantifies  the  nonlinear  and  hysteretic  behavior  of  the 
materials  through  a  broad  temperature  range. 

The  physical  mechanisms  which  provide  relaxor  ferroelectric  materials  with  their  characteristic  behavior 
and  a  general  form  of  the  constitutive  equations  which  quantify  this  paper  are  summarized  in  the  remainder 
of  this  section.  In  Section  2,  the  isothermal  anhysteretic  and  hysteresis  models  are  summarized  to  illustrate 
the  framework  of  the  model.  The  fully  temperature-dependent  model  is  then  developed  in  Section  3 
and  illustrated  through  fits  to  PMN  and  PMN-PT-BT  data  in  Section  4.  These  examples  demonstrate 
the  ability  of  the  model  to  quantify  both  the  hysteretic  and  anhysteretic  behavior  of  the  materials  for 
temperatures  ranging  from  133°  K  to  313°  K. 

1.1  Physical  Mechanisms  for  Relaxor  Behavior 

Normal  ferroelectric  materials,  including  the  piezoelectric  compound  PZT,  exhibit  a  sharp  ferroelectric 
to  paraelectric  phase  transition  and  are  operated  at  temperatures  well  below  this  transition.  In  contrast, 
relaxor  ferroelectrics  exhibit  a  diffuse  phase  transition  over  a  broad  temperature  range  [5,  7,  23].  Since 
relaxors  are  typically  employed  at  temperatures  within  this  transition  region,  the  material’s  dielectric 
behavior  strongly  depends  on  the  temperature  and  frequency  of  the  applied  electric  load  as  illustrated  in 
Figure  2.  The  transition  temperature  is  often  tailored  to  a  specific  application  by  mixing  PMN  with  lead 
titanate  (PT)  and  additional  minor  dopants,  such  as  barium  titanate  (BT)  and  strontium  titanate  (ST). 

The  temperature-dependence  of  relaxors  is  manifested  in  both  the  degree  of  hysteresis  and  the  satu¬ 
ration  exhibited  by  the  materials  as  depicted  in  Figure  1.  To  be  applicable  throughout  the  temperature 
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Figure  2.  Dielectric  response  of  a  PMN-PT-BT  (10%  PT,  3%  BT)  at  frequencies  ranging  from  0.1  kHz 
to  100  Hz. 

range,  models  must  accommodate  both  of  these  phenomena.  Furthermore,  Viehland  et  al.  [25,  26]  have 
demonstrated  that  relaxors  deviate  from  the  Curie- Weiss  law  at  temperatures  above  the  temperature  of 
dielectric  maximum,  Tm,  and  that  they  exhibit  non-Debye  relaxation  cumulating  in  polarization  freezing 
at  a  temperature  below  Tm. 

To  motivate  the  underlying  mechanisms,  it  is  necessary  to  consider  aspects  of  the  stoichiometry  of 
these  materials.  The  PMN  compounds  have  the  ABO3  shown  in  Figure  3.  The  A  sites  are  lead  anions, 
while  the  B  sites  are  a  mixture  of  low  valence  cations  and  high  valence  cations  (1/3  Mg+2  and  2/3  Nb+5 
in  the  case  of  PMN).  The  spatial  ordering  of  the  multiple  B-site  cations  leads  to  the  diffuse  transition 
of  a  relaxor.  Smolensky  postulated  a  dispersion  of  cations  which  creates  Mg:Nb  rich  ferroelectric  regions 
surrounded  by  pure  Nb  sites  [23].  Each  region  has  a  local  phase  transition,  so  the  aggregate  possesses 
a  spectrum  of  phase  transitions.  Cross  [5]  later  proposed  the  more  complete  picture  of  relaxor  behavior 
depicted  in  Figure  4.  He  postulated  that  nanoscale  size  (3-5  nm)  polar  regions  form  with  decreasing 
temperature.  These  small  regions  have  a  low  energy  barrier  so  they  fluctuate  with  thermal  agitation  in 
a  manner  analogous  to  superparamagnetism.  Each  region’s  dipole  moment  has  random  orientation  so 
that  the  macroscopic  polarization  is  zero  but  the  system’s  RMS  dipole  moment  is  nonzero.  Initially,  the 
relaxor  obeys  the  Curie- Weiss  relationship  with  a  very  high  Curie  temperature  (398°  K  for  PMN).  With 


Figure  3.  Atomic  structure  of  a  perovskite  relaxor  ferroelectric  A(B’,B”)03-  The  B’  and  B”  sites  indicate 
the  cation  ordering  for  the  “random  layer”  model. 
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Increasing  Temperature 


Ferroelectric:  The  micropo- 
lar  clusters  interact  and  ma¬ 
terial  exhibits  remanent  polar¬ 
ization. 


Superparaelectric:  Microp- 
olar  regions  develop  but  have 
a  random  orientation  due  to 
thermal  excitation. 


Paraelectric:  Micro-regions 
are  nonpolar  and  material  ex¬ 
hibits  zero  macroscopic  polar¬ 
ization. 


Figure  4.  Polar  mechanisms  in  relaxor  ferroelectric  materials  as  postulated  by  Cross  [5].  In  the  super¬ 
paraelectric  phase,  the  exhibits  a  nonzero  RMS  dipole  moment  but  has  zero  macroscopic  polarization. 


further  cooling,  the  thermal  fluctuations  decrease  and  interaction  between  regions  freezes  the  polarization 
in  a  manner  similar  to  magnetic  spin  glasses.  The  material  deviates  from  the  Curie- Weiss  relationship 
and  possesses  a  strong  frequency  dependence.  Eventually,  the  material  develops  a  permanent  remanent 
polarization  and  large  domains  develop  as  in  a  normal  ferroelectric. 

The  ordering  of  the  B-site  cations  plays  a  crucial  role  in  defining  the  diffuse  nature  of  the  transition  and 
two  cation  ordering  models  have  been  proposed.  Diffraction  studies  indicate  that  the  polar  regions  have  a 
1-to-l  distribution  of  two  different  cation  sites  arranged  in  a  face  center  lattice.  These  two  sites  are  denoted 
as  B’  and  B”  as  depicted  in  Figure  3.  In  the  “space  charge”  model  for  PMN,  it  is  assumed  that  ordered 
regions  consisting  of  Mg  exclusively  occupy  the  B’  sites  while  Nb  exclusively  occupies  the  B”  sites  [8,  26]. 
The  ordered  regions  thus  carry  a  net  negative  charge  which  should  inhibit  domain  growth.  In  the  “random 
layer”  model  [1,  3],  Mg  and  Nb  in  a  2-to-l  ratio  randomly  populate  the  B’  sites,  while  Nb  exclusively 
occupies  the  B”  sites.  In  this  model,  the  ordered  regions  are  charge-balanced  so  domains  should  coarsen 
during  heat  treatment.  Recent  experiments  by  Akbas  and  Davies  [1]  have  shown  that  domain  growth 
does  occur  during  annealing  in  Pb(Mgj/3  Ta2/3)03  (PMT),  which  supports  the  “random-layer”  model  of 
the  B-site  cation  distribution.  It  is  this  latter  model  which  we  employ  in  Section  3  when  quantifying  the 
macroscopic  behavior  of  the  materials. 

1.2  General  Constitutive  Relations 

We  summarize  here  the  general  form  of  the  constitutive  relations  for  relaxor  ferroelectric  materials.  The 
relation  for  the  direct  electrostrictive  effect  is  posed  in  terms  of  a  general  function  which  is  then  specified  in 
subsequent  sections.  To  simplify  the  discussion,  we  focus  on  the  1-D  case.  This  provides  models  which  are 
adequate  for  a  number  of  applications  utilizing  relaxors  and  provides  a  framework  which  can  be  extended 
to  multiaxial  models  if  the  application  warrants. 

To  define  an  appropriate  Gibbs  function  for  relaxor  materials,  we  assume  that  the  material  is  linearly 
elastic,  is  electromechanically  coupled  through  electrostrictive  effects,  and  that  it  exhibits  a  nonlinear 
and  potentially  hysteretic  relation  between  the  input  field  E.  the  macroscopic  polarization  P,  and  the 
operating  temperature  T.  For  polycrystalline  materials  which  are  approximately  isotropic,  an  appropriate 
Gibbs  function  is 

G  =  -±- a 2  -  Q33P2a  +  Q0{P,T) 
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where  a  denotes  the  axial  stress,  Y  is  the  open  circuit  elastic  modulus,  Q 33  is  the  longitudinal  electrostrictive 
coefficient  and  Qq  is  the  stress-free  Gibbs  free  energy.  The  use  of  the  thermodynamic  relations 


e  = 


dG 

da 


then  yields  the  general  constitutive  relations 


e  =  y°  +  Q33-P2 
E  =  —  2Q33P a  +  P(P,  T ) 


(1) 


specifying  the  strain  e  and  field  E  as  functions  of  the  stress,  polarization  and  temperature.  The  nonlinear 
map  J~ .  which  quantifies  the  stress-free  dielectric  response  of  the  material,  will  be  specified  in  subsequent 
sections. 


2  Isothermal  Hysteresis  Model 

To  motivate  the  strategy  used  to  construct  the  temperature-dependent  hysteresis  model  and  constitutive 
relations,  we  summarize  first  the  model  developed  in  [19,  20]  for  fixed  temperature  applications.  The  fixed 
temperature  model  is  constructed  in  two  steps:  (1)  the  characterization  of  the  hysteresis-free  (anhysteretic) 
relation  between  the  input  field  E  and  polarization  P,  and  (2)  the  modeling  of  hysteresis  losses  through 
the  quantification  of  energy  required  to  bend  and  translate  domain  walls  pinned  at  inclusions  in  the 
material.  We  note  that  while  the  anhysteretic  model  can  be  multivalued  for  certain  parameter  choices,  the 
transition  between  saturation  values  is  typically  steeper  than  that  observed  in  polycrystalline  materials 
as  temperatures  approach  the  Curie  point.  The  inclusion  of  domain  wall  loss  mechanisms  modulates  the 
transition  and  provides  accurate  predictions  for  the  polarization  at  a  range  of  fixed  temperatures.  An 
important  attribute  of  the  model  is  the  property  that  parameters  can  be  asymptotically  correlated  to 
physical  properties  of  the  data  including  the  coercive  field  and  remanence  point,  and  pertinent  asymptotic 
relations  from  [21]  are  summarized  at  the  end  of  the  section  to  motivated  components  of  the  temperature- 
dependent  model  developed  in  Section  3. 

Anhysteretic  Polarization 

The  anhysteretic  polarization  at  a  given  field  level  E  is  defined  to  be  the  polarization  that  results 
when  a  decaying  AC  field  centered  at  E  is  applied  to  the  material.  For  a  fixed  stress,  it  represents  an 
equilibrium  between  the  thermal  and  electrostatic  energies  and  can  thus  be  modeled  using  Boltzmann 
principles.  Under  the  assumption  that  dipoles  can  orient  only  in  the  direction  of  the  applied  field,  or 
opposite  to  it,  the  balance  of  energy  yields  the  Ising  spin  model 

Pan  =  PoAAtanh  (2) 

where  po  is  the  dipole  moment  of  a  single  cell  and  M  is  the  total  number  of  moments  per  unit  volume 
{M  =  N/V).  The  thermal  energy  is  quantified  by  kisT  where  ku  is  Boltzmann’s  constant  and  T  is  the 
temperature. 

To  complete  the  model  for  the  anhysteretic  polarization,  it  is  necessary  to  incorporate  interactions  due 
to  neighboring  dipoles.  These  effects,  which  are  analogous  to  the  Weiss  fields  in  magnetic  materials,  are 
modeled  under  the  assumption  that  the  effective  field  at  the  domain  level  has  the  form 

Ee  =  E  +  otP  +  Pa  .  (3) 
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The  component  aP  incorporates  the  effects  of  interdomain  coupling  while  Pa  quantifies  the  electrome¬ 
chanical  coupling  due  to  applied  stresses.  Finally,  to  guarantee  that  the  model  saturates  to  the  saturation 
polarization  Ps,  it  is  necessary  to  define 

Ps  -  poAf  ■ 

For  regimes  in  which  the  electromechanical  coupling  is  adequately  modeled  by  the  bilinear  relation  Pa  — 
‘IQx^Pa.  the  field  can  then  be  expressed  as 


E  —  —2Q^Pa 


aP  +  -TT-pp- arctanh  (P/Ps ) . 
y  -is 


The  function  T  in  (1)  is  this  given  by 


P(P.T)  = 


—aP  + 


kBNT 

VPS 


arctanh(P/Ps) . 


(4) 

(5) 


For  implementation  purposes,  it  is  advantageous  to  express  the  anhysteretic  polarization  as 


Pan  -  Pstanh 


(6) 


where  a  —  like  a  and  Ps.  is  a  macroscopic  parameter  which  is  estimated  using  either  the  asymptotic 
relations  developed  in  [21]  or  through  a  least  squares  fit  to  data.  We  note  that  while  a  form  of  temperature- 
dependence  is  incorporated  in  the  parameter  a,  it  does  not  provide  a  complete  characterization  of  the 
thermal  behavior  of  relaxor  ferroelectrics  through  their  full  operating  range.  The  parameters  a  and  Ps 
also  vary  with  temperature,  and  appropriate  expressions  quantifying  the  thermal  dependence  of  all  three 
parameters  are  developed  in  Section  3. 

For  certain  materials,  the  assumption  that  dipoles  orient  either  with  or  against  the  field  is  too  stringent 
and  produces  models  which  saturate  too  quickly.  The  less  restrictive  assumption  that  dipoles  can  orient 
uniformly  yields  the  Langevin  relation 


p  —  p 
J  an  —  s 


(7) 


when  the  thermal  and  electrostatic  energies  are  balanced  through  Boltzmann  relations.  The  anhysteretic 
model  (7)  can  provide  slightly  improved  accuracy  when  characterizing  certain  materials  but  it  cannot  be 
inverted  to  yield  polarization-based  constitutive  relations  comparable  to  (4).  Both  the  Ising  spin  expression 
(6)  and  Langevin  relation  (7)  have  been  employed  when  modeling  the  anhysteretic  behavior  of  ferroelectric 
and  ferromagnetic  materials. 


Hysteresis  Model 

The  second  component  of  the  model  quantifies  the  energy  required  to  translate  and  bend  domain  walls 
pinned  at  inclusions  in  the  material.  As  detailed  in  [19,  20],  the  irreversible  changes  in  polarization  due  to 
domain  wall  translation  are  incorporated  through  the  quantification  of  the  energy 

£pin(P )  =  k  [  dP 

Jo 

required  to  reorient  dipole  moments.  The  parameter  k  is  defined  by  k  —  where  n  denotes  the  average 
density  of  pinning  sites,  p  is  an  average  dipole  moment,  and  (£w)  is  the  average  energy  for  180°  domain 
walls.  While  k  was  considered  constant  in  the  fixed-temperature  analysis  of  [19,  20,  21],  it  is  clearly 
temperature-dependent  due  to  the  dependence  of  {£v)  on  the  thermal  energy. 
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(8) 


A  balance  of  electrostatic  energy  then  yields  the  relation 

(]P  ~  P  —  P 

irr  g  1  an  1  irr 

dE  kS  OL  ( Pan  Pirr ) 

for  the  irreversible  polarization.  The  parameter  8  =  sign  (A)  ensure  that  the  energy  required  to  break 
pinning  sites  always  opposes  changes  in  the  polarization.  The  physical  observation  that  polarization  changes 
after  a  reversal  in  field  direction  are  reversible  motivates  the  incorporation  of  the  parameter  <5  which  is 
defined  to  1  if  {dE  <  0  and  P  >  Pan}  or  {dE  >  0  and  P  <  Pan}  and  is  0  otherwise.  The  solution  of  (8) 
quantifies  the  irreversible  changes  in  polarization  due  to  domain  wall  translations. 

The  second  component  of  the  polarization  is  the  reversible  polarization  which  models  the  effects  of 
domain  wall  bending.  To  a  first  approximation,  this  is  modeled  by  the  relation 


Prev  —  ('{Pan  Pirr ) 


(9) 


where  c  is  a  parameter  which  is  estimated  for  a  given  material  [19,  20,  21]. 

The  total  polarization  is  given  by 

P  =  Pirr  +  Prev  ■  (10) 

To  implement  the  model,  the  effective  field  for  a  given  field  and  irreversible  polarization  level  is  computed 
using  (3).  This  effective  field  value  is  then  employed  in  either  (6)  or  (7)  to  compute  the  corresponding 
anhysteretic  polarization.  The  subsequent  irreversible  polarization  is  determined  by  numerically  integrating 
(8).  The  total  polarization  is  then  specified  by  (10). 

To  specify  the  functional  T  in  the  general  constitutive  relations  (1),  it  is  necessary  for  the  field  to 
be  expressed  as  a  function  of  the  polarization.  This  can  be  accomplished  only  through  the  solution  of  a 
differential  or  integral  equation  since  this  is  the  form  of  the  model  used  to  quantify  the  evolving  hysteresis. 
The  construction  of  a  consolidated  differential  equation  which  facilitates  inversion  is  detailed  in  [20]  and 
details  regarding  the  solution  of  the  inverted  problem  is  described  in  [17,  18,  22]  in  the  context  of  the 
analogous  magnetization  model.  Hence  while  one  cannot  provide  an  algebraic  constitutive  relation  for  the 
polarization-based  field  for  this  hysteresis  model,  the  map  between  P  and  E  is  well-defined  as  the  solution 
to  an  evolution  equation.  The  determination  of  the  inverse  through  the  solution  of  a  complementary 
differential  equation  also  facilitates  the  construction  of  inverse  compensators  for  linear  control  design. 


Asymptotic  Parameter  Relations 

The  previous  model  requires  the  parameters  a,  a,  c,  Ps  and  k.  Algorithms  for  determining  initial  es¬ 
timates  of  these  parameters  using  measured  values  of  the  remanence  polarization,  coercive  field,  peak 
polarization  and  field,  and  the  differential  susceptibilities  at  these  points  were  developed  in  [21].  We  sum¬ 
marize  here  those  asymptotic  relations  which  pertain  to  the  development  of  the  temperature-dependent 
model. 

Let  Xan  and  Xc  respectively  denote  the  differential  susceptibilities  of  the  anhysteretic  polarization  at 
E  =  0  and  the  total  polarization  at  E  =  Ec  as  depicted  in  Figure  5.  The  slopes  of  the  polarization  curve 
before  and  after  tip  reversal  are  denoted  by  and  xm-  As  developed  in  [21],  estimates  of  the  parameters 
a,  c  and  k  are  given  by  the  asymptotic  relations 


Ps 

a  3  Xan 

(11a) 

r  —  Am 

°  —  + 

Xm 

(lib) 

k  ra  Ec  (For  soft  materials). 

(11c) 
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Figure  5.  Hysteresis  curve  with  differential  susceptibilities  employed  for  parameter  determination. 

The  saturation  polarization  Ps  can  quite  accurately  be  estimated  directly  from  high  drive  level  data  while 
the  slope  Xan  is  approximated  by  the  slope  at  coercivity  Xan  ~  Xc-  Details  concerning  these  relations  as 
well  as  the  complete  algorithm  for  estimating  the  full  parameter  set  are  provided  in  [21]. 

3  Temperature-Dependent  Hysteresis  Model 

From  a  physical  perspective,  the  parameters  a,  a,  Ps,  k  and  c  are  temperature-dependent  so  the  model  (10) 
is  valid  only  in  operating  regimes  in  which  temperatures  exhibit  minimal  variation.  In  this  section  we 
address  the  temperature-dependence  through  the  consideration  of  the  relaxor  ferroelectric  as  an  aggregate 
of  micropolar  regions  having  differing  densities  of  Mg  cations  and  hence  differing  ferroelectric  properties. 
The  development  of  the  anhysteretic  relations  is  motivated  by  Glazounov  et  al.  [7]  and  follows  the  theory 
developed  by  Horn  and  Shankar  [13],  The  temperature-dependence  in  the  pinning  losses  is  incorporated 
through  energy  arguments  analogous  to  those  employed  in  the  theory  of  dislocation  plasticity  [14].  The 
combined  model  can  be  used  to  characterize  hysteresis  in  relaxor  ferroelectrics  throughout  a  wide  range  of 
temperatures  encountered  in  operation. 

3.1  Anhysteretic  Polarization 

As  noted  previously,  relaxor  ferroelectrics  are  assumed  to  consist  of  micropolar  regions  whose  transition 
temperatures  and  ferroelectric  properties  are  determined  by  variations  in  the  concentration  of  cations 
(Mg2+  and  Nb°+  in  PMN)  which  comprise  the  material.  The  anhysteretic  component  of  the  model  is 
constructed  in  two  steps.  In  the  first,  we  consider  the  thermodynamics  of  the  micropolar  regions  in  order 
to  predict  the  saturation  polarization  and  density  of  Mg  cations  as  a  function  of  temperature.  Secondly,  the 
interaction  between  regimes  is  quantified  to  construct  a  macroscopic  model  which  predicts  the  dielectric 
properties  and  polarization  of  the  material  for  temperatures  below  the  freezing  temperature. 

For  modeling  purposes,  we  make  the  following  assumptions  which  are  depicted  in  Figure  6. 

(i)  Each  region  contains  a  total  of  Nq  Mg  cations. 

(ii)  The  density  of  Mg  cations  X  is  uniform  throughout  each  region  but  varies  between  regions.  As  a 
result  of  Assumption  (i),  the  volume  V%  is  each  region  is  inversely  proportional  to  X,t. 
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Figure  6.  Relaxor  ferroelectric  consisting  of  an  aggregate  of  micropolar  regions  which  satisfy  Assumptions 
(i)  and  (ii).  In  this  depiction,  there  are  Nr  —  4  regions  and  a  total  of  N  —  NrNq  Mg  cations  in  the  material. 


(iii)  The  strength  of  interactions  between  regions  is  directly  proportional  to  the  density  of  nearest  neigh¬ 
bors. 

(iv)  Below  its  local  transition  temperature,  each  region  acts  as  a  single  dipole  and  exhibits  a  remanent 
polarization. 


3.1.1  Thermodynamics  of  Micropolar  Regions 

We  consider  first  the  behavior  of  the  micropolar  regions  which  is  determined  by  the  density  of  Mg  cations 
which  occupy  B’  sites.  Under  the  assumption  that  the  regions  consists  of  a  lattice  of  cells  as  depicted  in 
Figure  3,  the  dipole  moment  p  for  the  lattice  is  given  by 

N0 

p  =  po'52si 

i—  1 

where  po  is  the  dipole  moment  for  a  single  cell  and  Nq  denotes  the  number  of  Mg  cations  in  the  region. 
For  PMN,  only  cells  containing  a  Mg  cation  possess  an  associated  spin  s*  =  ±1  where  the  sign  depends 
upon  the  cell’s  orientation.  A  correlation  between  cells  is  provided  by  the  average 


R  = 


N0 


or  equivalently,  R  =  7^. 

As  noted  in  Assumptions  (i)  and  (ii),  we  assume  that  the  regions  all  have  the  same  number  Nq  of 
Mg  cations  but  differing  densities  and  hence  volumes.  Within  a  region,  we  let  N+  and  iV_  respectively 
denote  the  number  of  positive  and  negative  cells.  By  noting  that  Nq  —  N+  +  JV_,  the  cell  numbers  can  be 
expressed  in  terms  of  the  correlation  through  the  relations 


N+  =  ^-(l  +  R) 

W_  =  ^°(l-R). 


(12) 


9 


The  density  X  is  defined  as  the  number  of  neighboring  B’  sites  that  contain  a  Mg  cation.  As  indicated 
by  Figure  3,  X  can  range  from  1  to  12.  Furthermore  X  is  assumed  to  be  uniform  throughout  the  region 
in  accordance  with  Assumption  (ii). 

As  will  be  detailed  later,  the  microregions  have  a  distribution  of  Curie  temperatures  due  to  the  differing 
cation  densities  in  the  regions.  To  determine  these  densities,  we  note  that  B’  sites  have  a  ratio  of  2/3  Mg 
cations  to  1/3  Nb  ions.  For  a  Mg  populated  B’  site,  each  of  the  12  adjacent  B’  sites  has  a  65.4%  probability 
of  being  populated  by  a  Mg  cation  and  34.6%  probability  of  being  populated  by  Nb.  The  probability  that 
the  site  is  surrounded  by  X  Mg  cations  is  computed  by  assuming  a  binomial  distribution 

1 21 

B(X)  =  xl{uix)l( 0-654)x(0.346)12-x  . 

To  obtain  a  continuous  distribution  appropriate  for  subsequent  integration,  this  was  then  fit  with  the 
normal  distribution 

IIpO  =  — ^e-(x-a)2/262  (13) 

by/  2ir 

having  a  mean  a  —  7.85  and  standard  deviation  b  —  1.65.  Finally,  the  distribution  is  normalized  to  the 
interval  [0, 1]  through  the  introduction  of  the  variables 

c  =  -  C  =  - 
s  12  ’  s  12 

to  obtain 

n(C)  =  ^=e-144^-«2/262  .  (15) 

The  distributions  (13)  and  (15)  will  be  employed  in  subsequent  discussion  to  compute  ferroelectric  prop¬ 
erties  of  the  material  which  depend  upon  cation  densities  in  the  micropolar  regions. 

To  compute  the  internal  energy  due  to  spin  conversion,  we  consider  first  the  energy  required  to  convert 
an  individual  spin  given  a  density  X  of  neighboring  Mg  cations.  Letting  4>o  denote  the  energy  required  to 
convert  the  spin  at  a  single  site  when  the  system  is  completely  ordered  ( R  —  ±1),  the  energy  required  to 
convert  a  positive  spin  to  negative  (and  conversely)  is  given  by 


(16) 


We  point  out  that  these  energy  expressions  are  derived  under  the  assumption  that  dipoles  interact  only 
with  adjacent  neighbors. 

The  change  in  the  internal  energy  U  due  to  spin  conversion  can  then  be  expressed  as 

dU  =  .  (17) 

By  utilizing  the  energy  relations  (16)  and  noting  from  (12)  that  dN+  —  4/dR,  —  —dN_.  the  expression 
(17)  can  be  integrated  to  obtain  the  relation 

U  =  ^M(i  -  R?)  +  U0r 

where  Uq  denotes  the  energy  for  the  completely  ordered  state,  for  the  internal  energy. 


X  N, 


=  TT7  " 


12  N0 
X  N- 


-<f>n 


12  N0 


o  • 


(14) 
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For  stress  free  conditions,  the  total  Gibbs  free  energy  T  for  the  microregion  is  related  to  the  internal 
energy  by 

T  =  U-  STV0 

where  S  denotes  the  entropy  and  Vo  is  the  volume  of  the  microregion.  From  statistical  mechanics,  the 
entropy  per  unit  volume  for  the  system  can  be  expressed  as 


S 


kBN o 

TT 


ln[io(i?)] 


UbNq 

TT 


In 


IY  N°  V 

kBN0  1 

[  JVo!  1 

K  w+  \ 

“  Vo  ln 

N+\N_\_ 

where  w(R)  denotes  the  number  of  possible  spin  combinations  that  yield  the  correlation  R.  The  entropy 
expression  can  be  simplified  to 


S  = 


kbn0 

Vo 


ln(2)  -  Ml  +  R)~  hJ*  Ml  -  R) 


which  yields  the  free  energy 

d'o^VoC 


r  = 


(1  -R2)  + 


kBNoT 


R  In 


1  +  R 
1  -R 


-ln(l  -i?2) 


+  Uq  +  kBN, o  ln(2) 


for  the  region.  We  note  that  the  expression  for  F  incorporates  the  ordering  of  neighboring  moments  reflected 
in  the  first  term  and  the  tendency  for  disorder  due  to  thermal  agitation  as  reflected  in  the  second  term. 
The  electric  field  Eioc  in  the  polar  region  is  determined  by 

9r  r 

p0N0Eioc  =  —  = - j - h  TA;Barctanh(i?) . 

In  the  absence  of  a  local  field,  the  equation  can  be  expressed  as 


R  = 


-arctanh(f?) 


(18) 


where 


T  — 
±cur  — 


C£o 
2 kB 


(19) 


denotes  the  local  Curie  temperature  for  that  region.  To  motivate  the  definition  (19)  and  this  terminology, 
it  is  noted  that  (18)  has  a  single  solution  for  T  >  Tcur  and  three  solutions  comprised  of  two  stable  and  one 
unstable  solution  for  temperatures  T  <  Tcur.  Hence  a  pitchfork  bifurcation  occurs  at  T  =  Tcur  which  is 
consistent  with  a  second-order  phase  transition. 

The  expression  (19)  quantifies  the  Curie  temperature  for  a  single  micropolar  region.  The  macroscopic 
Curie  temperature  for  the  material  is  specified  as  the  average  Curie  temperature  for  the  regions  or  the 
Curie  temperature  of  the  mean, 

C$o 


T  = 

cur 


2  kB 


(20) 


The  behavior  of  each  polar  regions  can  then  be  expressed  as 


$0 

PqNqEioc  =  — 


CT 

— C R  +  = — arctanh(i?) 

T  cut 


(21) 


in  the  absence  of  interactions  between  regions. 
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3.1.2  Macroscopic  Relations 


The  previous  discussion  focuses  on  the  thermodynamics  within  a  micropolar  region  and  does  not  address 
interactions  between  regions.  As  detailed  in  [20],  however,  such  interactions  can  significantly  affect  the  ef¬ 
fective  field  present  within  the  material  and  we  focus  here  on  the  development  of  commensurate  constitutive 
relations  for  the  aggregate  material  which  incorporates  these  effects. 

The  constitutive  relation  (4)  quantifies  the  dependence  of  the  field  on  the  polarization  for  fixed  temper¬ 
atures.  To  incorporate  thermal  dependence,  it  is  necessary  to  quantify  the  effects  of  temperature  on  the 
number  of  dipoles  N  and  the  coupling  coefficient  a.  As  depicted  in  Figure  6,  N  is  taken  to  be  N  —  NRNo 
where  NR  and  No  respectively  denote  the  number  of  micropolar  regions  and  the  number  of  Mg  cations  in 
each  region.  To  quantify  a,  we  employ  Assumption  (iii)  from  Section  3.1  and  assume  that  the  strength 
of  interactions  is  proportional  to  the  number  of  neighbors  Nd  within  a  distance  d  of  a  cluster.  Under 
stress-free  anhysteretic  conditions,  we  then  propose  the  constitutive  relation 


E  = 


®NdNR  kBTNRNo 


2VP} 


-P  + 


VP* 


arctanh  (  — 


(22) 


for  relaxor  materials.  Here  V  denotes  the  volume  of  material  while  $  is  a  relative  measure  of  the  interaction 
energy  between  polar  regions.  We  note  that  the  second  term  on  the  right  side  of  (22)  is  proportional  to 
the  function  employed  by  Glazounov  et  al.  [7]  for  non-interacting  clusters. 

Under  the  assumption  that  the  clusters  are  uniformly  distributed,  the  number  within  a  distance  d  can 
be  expressed  as 

M  4  ird3NR 

X d  — 


3U 


With  the  definition 


Tf  = 


_  2ir$d3NR(Tf) 


3  kBV 


for  the  freezing  temperature,  the  temperature-dependent  constitutive  relation  between  the  field  and  polar¬ 
ization  can  then  be  expressed  as 


E(T)  = 


2-K<&d?NR 

3V 


-Nr 

VP'i 


P  + 


NoNR(Tf)  T 

PsV  2> 


arctanh 


(23) 


We  now  employ  the  theory  developed  in  Section  3.1.1  for  the  thermodynamics  of  the  micropolar  regions 
to  express  Ps  and  NR  in  terms  of  the  operating  temperature  and  parameters  which  must  be  estimated  for 
a  specific  material. 

The  saturation  polarization  occurs  when  the  dipole  moments  in  all  the  microregions  are  aligned  with 
the  field.  The  degree  of  alignment  for  microregions  is  quantified  by  the  correlation  R(T,X)  specified  by 
(21)  while  the  magnitude  is  dependent  upon  the  number  n(X)  of  Mg  cations  having  X  nearest  neighbors. 
Recalling  that  po  denotes  the  dipole  moment  of  a  single  cell,  the  saturation  polarization  can  be  expressed 
as 

Ps(T)  =  ^  I'2  \R(T,X)\p0n(X)dX  . 

V  Jo 

Since  N  denotes  the  total  number  of  Mg  cations  in  the  volume  V.  we  can  express  n(X)  —  NII(X). 
Normalization  to  the  unit  interval  then  yields  the  expression 

ps(T)  =  Ps  /VcGOintcR-  (24) 

Jo 


The  distribution  n  is  computed  using  (15)  and  the  parameter  Ps  —  is  estimated  either  directly  from 
the  data  or  through  a  least  squares  fit  to  the  data.  The  correlation  R  is  computed  through  iterative  solution 


12 


of  the  relation 


(3Ei0C  —  —(R  + 


-arctanh(P) 


where  (3  =  2p^°  denotes  the  average  degree  to  which  local  fields  turn  the  polarization  and  hence  change 
the  correlation;  see  (21).  We  note  that  at  low  drive  levels,  this  contribution  is  small  and  one  can  employ 
(3  =  0.  For  high  drive  levels,  however,  Ps  is  more  accurately  approximated  by  treating  (3  as  a  parameter 
to  be  estimated  and  solving  (25)  with  Eioc  =  Emax. 

The  number  of  micropolar  regions  Nr,  at  a  fixed  temperature  T,  is  computed  by  determining  those 
regions  with  nonzero  correlation  R,  or  from  (25),  those  regions  in  which  T  <  Tcur .  Since  the  number  of 
regions  with  density  X  is  given  by  n{X)/No  =  NU(X)  /No,  the  total  number  of  regions  in  the  volume  is 

Nr{T)  _  N  P  _  12u(QdC  • 

J  CT /Tcur 


V  VNo  Jct/Tcut  v"'  " 

The  use  of  (24)  and  (26)  in  (23)  yields  the  temperature-dependent  constitutive  relation 
„  -2vr $d3  /  N  \2  I?(T )  „  2t r$d3W0  (  N  \2  U  P  (T)  ,  (  P 


I2(T)  2-K^d3  Np 

P?(T)  3 


2  h_  M^Lrctanh 
TfPs{T ) 


Ps(T) 


where 


h  (T)  =  f! 

■k'i 

w:  _ 

kif/Tc 


i2n(C)d( 


1211(C)  d( 


ps(t)  =  ps  fnmcmodc. 

Jo 


For  implementation  purposes,  we  formulate  (27)  as 

E  =  a^prP  +  a— ^^-arctanh  (— ^—)  (29) 

PS2(T)  TfPs{T)  \PS(T)  J  y  1 

where  a  and  a  are  parameters  which  are  estimated  for  a  given  material  through  a  least  squares  fit  to  data. 
For  the  subsequent  hysteresis  model,  the  constitutive  relation  is  formulated  as 


E  =  a(T)P  +  a(T)arctanh 


Ps(T) 


with  the  temperature-dependent  parameters 

a(T)  =  a 


km 

P2(T) 


m  _~h  h(T) 
(  }  TfPs(T) 


PS(T)  =  PS  f'ulRiTXmOdC 

Jo 


where  Ps  =  must  be  estimated  for  a  given  material.  The  resulting  model  for  the  anhysteretic  polar¬ 
ization  is 


Pan{T)  =  PS{T)  tanh 


E  +  a(T)Pan{T) 
a(T) 
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3.2  Domain  Wall  Losses 


The  anhysteretic  model  (31)  quantifies  the  polarization  which  would  occur  in  the  absence  of  domain  wall 
pinning.  For  certain  parameter  choices,  this  model  produces  a  form  of  hysteresis  similar  to  that  observed  in 
a  single  crystal  which  undergoes  domain  switching.  However,  the  transition  from  remanence  to  the  coercive 
point  is  steeper  than  that  observed  in  most  operating  regimes.  Furthermore,  the  parameters  obtained  at 
temperatures  rising  to  the  material’s  freezing  point  typically  yields  a  single-valued  anhysteretic  curve 
whereas  hysteresis  is  still  present  in  the  E-P  relation.  These  phenomena  are  attributed  to  energy  losses 
due  to  domain  wall  pinning. 

As  summarized  in  Section  2,  the  quantification  of  the  irreversible  and  reversible  polarizations,  respec¬ 
tively,  requires  the  parameters  k  and  c.  We  focus  here  on  the  quantification  of  the  temperature-dependence 
in  k  since  it  directly  affects  the  decrease  in  coercivity  observed  as  temperatures  increase  to  the  freezing 
point.  The  parameter  c  is  taken  as  constant  since  its  variation  with  temperature  appears  to  be  secondary. 
If  further  refinement  is  required,  analysis  of  the  type  developed  for  k  can  be  extended  to  include  the 
reversible  effects  quantified  by  c. 

For  fixed  temperature,  it  is  noted  in  Section  2  that  through  the  definition  k  —  rt^  ■ ,  the  parameter  k 
depends  on  the  180°  domain  wall  energy  (£w)  required  to  move  between  potential  wells.  The  temperature 
dependence  in  k  is  attributed  to  the  temperature-dependence  of  the  potential  wells  and  hence  (£1).  To 
quantify  the  reduction  in  energy  required  to  move  between  wells  at  increasing  temperatures,  we  employ 
the  asymptotic  relation  k  &  Ec  noted  in  (11c)  and  focus  initially  on  the  temperature-dependence  of  the 
coercive  field  Ec. 

The  influence  of  temperature  on  the  coercive  field  behavior  is  predicted  using  energy  arguments  drawn 
from  the  theory  of  dislocation  plasticity  [14]  and  is  similar  to  the  approach  employed  by  Chen  and  Wang 
[4],  In  this  theory,  energy  barriers  AG  due  to  local  resistive  fields  E  are  considered  as  the  mechanism  which 
inhibits  domain  wall  movement  and  hence  prohibits  domain  wall  movement  and  hence  provides  potential 
wells.  Using  energy  arguments  analogous  to  those  for  slip  dislocations,  the  application  of  Boltzmann  theory 
yields  the  modified  Arrhenius  equation 


dP 

dt 


=  AP0( 


m 


Po  exp 


-AG\ 
kBT  ) 


(32) 


for  the  polarization  rate.  Here  po,  A <7,  h  and  APo  respectively  denote  the  frequency  at  which  the  domain 
wall  vibrates,  the  distance  traveled  by  the  wall,  the  distance  between  walls,  and  the  local  change  in 
polarization  due  to  domain  wall  movement. 

To  specify  the  energy  barrier  AG,  we  consider  the  phenomenological  model 


AG  =  G0 


(f)7 


which  is  analogous  to  that  employed  in  the  theory  of  dislocation  plasticity  [14] .  The  parameter  Go  desig¬ 
nates  the  total  free  energy  required  to  overcome  the  obstacle.  The  parameters  p  and  q  quantify  the  shape 
of  the  obstacle  as  depicted  in  Figure  7.  For  the  analysis  which  follows,  the  choice  q  —  1  with  general  p 
proves  sufficiently  general,  and  in  this  case 

AG  =  G0(P1/p  -eVp) 

where  Go  =  Go/E 1/p .  We  note  that  in  their  theory,  Chen  and  Wang  employed  p  —  q  —  1. 

For  fixed  frequencies  which  are  sufficiently  low  to  permit  quasistatic  approximations,  this  yields 


dPdE 
dE  dt 


=  Pi  exp 


-G0  ( E1/p  -  E1/P^j 
kPl' 
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Figure  7.  The  influence  of  p  and  q  on  the  shape  of  short  range  obstacles  which  are  sensitive  to  thermal 
activation  (after  [14]). 


where  p\  —  AP^Agg^/h.  With  the  definition  p2  —  it  follows  that 

El,r=TwM^)+sl,r- 

Consideration  of  the  coercive  field  E  —  Ec  then  yields 

EVp  =  In  \p2Xc(T)\  +  E (33) 
Co 

where  Xc{T)  denotes  the  differential  susceptibility  at  Ec  (See  Figure  5).  The  susceptibility  will  be  temper¬ 
ature-dependent  since 

Xc{T)  ~  Xan{T)  ~ 

as  indicated  by  the  asymptotic  relation  (lib)  of  Section  2.2.  The  contribution  ln[p2Xc(F)]  is  dominated 
by  T,  however,  so  to  a  first  approximation,  we  can  represent  the  temperature  dependence  of  the  coercive 
field  by 

Ec  =  k(Tf  -  Tf  (34) 

where  Tj  is  the  freezing  temperature  where  remanence  disappears  and  A:  is  a  material-dependent  parameter. 
We  note  that  for  p  —  1,  the  expression  (34)  yields  a  linear  relation  between  the  coercive  field  and  operating 
temperature  which  is  consistent  with  the  data  plotted  in  Figure  2  of  [4]  for  a  variety  of  materials.  At  low 
temperatures,  however,  one  can  observe  a  quadratic  dependence  as  illustrated  for  PMN  in  Example  2  of 
the  next  section. 

Finally,  from  the  asymptotic  relation  k  ~  Ec  (see  (11c))  for  soft  materials,  we  consider  the  relation 

k(T)  =  k(Tf-T)*  (35) 

for  the  pinning  constant  at  temperatures  T  <  Tf  where  hysteresis  is  present.  For  temperatures  T  >  Tj, 
the  material  exhibits  no  hysteresis  so  k(T)  =0  and  the  irreversible  relation 

dP 

PirriT )  =  Pan (T)  -  6s0k(T)-^(T) ,  (36) 

derived  in  [20]  and  used  to  obtain  (8),  reduces  to  the  anhysteretic  model  (31). 
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3.3  Temperature-Dependent  Hysteresis  Model 

The  full  hysteresis  model  is  constructed  in  the  manner  summarized  in  Section  2  with  the  temperature- 
dependent  coefficients,  developed  in  Sections  3.1  and  3.2,  employed  in  the  models  for  the  anhysteretic, 
irreversible  and  reversible  polarizations.  The  parameters,  along  with  the  complete  model  are  summarized 
in  Algorithm  1  below.  Details  regarding  the  implementation  of  this  algorithm  are  provided  in  [20,  21]. 
Examples  illustrating  the  performance  of  the  model  are  provided  in  Section  4. 

Algorithm  1: 

(A)  Determine  Temperature-Dependent  Parameters 

(1)  At  a  specified  temperature,  estimate  a,a,k,  (3,c,p  and  Ps  =  poN/V  using  the  asymptotic  rela¬ 
tions  from  [21]  or  a  least  squares  fit  to  data.  Note  that  for  small  drive  level  applications  and 
moderate  temperature  ranges,  one  can  often  employ  (3  =  0  and  p  =  1. 

(2)  Compute  the  temperature-dependent  parameters 

<*(T)  =  ,  PS(T )  =  Ps  1 12|P(T,  C)|  n(C)  d( 

(37) 

a{T)  =  zTf¥Jyj  ’  M t)  =  k(Tf-Tf 

where 

h(T)=  f!  _  1211(C)  dC 

J  (,T /T  cur 

h  =  C  _  1211(0  d( 

Cl f  /T cur 

n(0  =  -J—e-lutt-02/2h2  )  ^  =  a/i2  ;  a  =  7.85,  b  =  1.65 

bV  2-k 

CT 

fdEmax  =  -( R  +  -=-arctanh(i?) 


(B)  Temperature-Dependent  Polarization 

(1)  Construct  Anhysteretic  Polarization 

Pan(T)  =  PS(T)  tanh 

(2)  Irreversible  Polarization 

dPirr 


E  +  a(T)Pjrr(T) 


dE 


-(T)=6- 


a(T) 

Pan(T)  ~  Pirr{T) 


(3)  Reversible  Polarization 

(4)  Total  Polarization 


kS-a[Pan(T)-Pirr(T)] 

Prev{T)  =  c[Pan(T)  -  Pirr{T )] 


P(T) 

—  Pirr  (T)  + 

Prev  (T) 
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4  Model  Validation 


The  model  derived  in  Section  3  quantifies  the  polarization  in  relaxor  ferroelectric  materials  as  a  function 
of  both  the  input  field  E  and  the  temperature  T.  To  illustrate  the  performance  of  the  model,  we  consider 
both  PMN  and  PMN-PT-BT  data  over  a  wide  range  of  temperatures  both  above  and  below  the  freezing 
temperatures.  For  each  set  of  data,  a  single  set  of  parameters  was  employed  throughout  the  full  temperature 
range,  and  the  temperature  and  electric  field  provided  the  sole  inputs  to  the  model.  The  results  illustrate 
that  through  the  incorporation  of  associated  physics,  the  model  characterizes  a  wide  range  of  dielectric 
behavior. 


Example  1  (PMN-PT-BT): 

We  consider  first  the  capability  of  the  model  to  characterize  and  predict  the  dielectric  behavior  of 
PMN-PT-BT  for  input  fields  having  a  peak  value  of  1  MV /m  at  temperatures  ranging  from  T  —  263°  K  to 
T  —  313°  K.  At  the  lower  temperature,  the  material  exhibits  ferroelectric  behavior  and  significant  hysteresis 
as  shown  in  Figure  8.  As  the  temperature  increases  through  the  freezing  temperature  (approximately 
288°  K),  the  hysteresis  vanishes  and  the  material  exhibits  a  single- valued  superparaelectric  relation  between 
the  applied  field  and  measured  polarization.  The  data  was  collected  from  a  stress-free  sample  (<r  =  0)  and 
the  input  field  was  cycled  at  1  Hz  to  maintain  quasi-static  conditions. 

The  parameter  values  a  —  1.5  x  107,a  =  7.65  x  107,  A:  =  4.0  x  103,PS  =  0.53,  (3  —  7  x  10~9,c  =  .3 
and  p  —  1  were  obtained  through  a  least  squares  fit  to  the  data.  The  freezing  temperature  was  specified 


-1  -0.5  0  0.5  1 


-1  -0.5  0  0.5  1 


-1  -0.5  0  0.5  1 


-1  -0.5  0  0.5  1 


Figure  8.  PMN-PT-BT  data  ( - )  and  the  domain  wall  model  ( - )  for  the  temperature  range 

T  —  263°K  to  T  —  313°K.  Abscissas:  electric  field  (MV/m),  ordinates:  polarization  (C/m2). 
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to  be  Tf  =  288°  K  and  the  Curie  temperature  was  taken  to  be  Tcur  =  313°  K.  We  reiterate  that  the 
temperature-dependence  in  the  parameters  a(T) ,  a(T) ,  k(T)  and  PS(T )  is  incorporated  solely  through  the 
relations  (37)  and  the  same  parameter  set  was  employed  through  the  full  range  of  temperatures  and  field 
inputs. 

The  model  fit  obtained  with  these  parameters  is  compared  with  the  measured  data  in  Figure  8.  It 
is  observed  that  the  model  very  accurately  quantifies  the  initial  hysteresis  at  T  =  263°  K,  the  decreased 
hysteresis  at  T  =  283°  K,  and  the  anhysteretic  behavior  at  higher  temperatures.  It  also  quantifies  the 
decrease  in  the  maximum  measured  polarization  which  occurs  as  temperatures  increase.  We  point  out  the 
model  parameters  were  chosen  to  optimize  the  model  predictions  throughout  the  full  temperature  range 
and  better  fits  can  be  obtained  for  specific  temperature  intervals  if  desired.  For  example,  the  parameter 
choices  5  =  1.4  x  107,  a  =  5.74  x  107,  k  =  5.0  x  103,  Ps  =  0.53  and  c  =  .5  yield  a  highly  accurate  model  fit 
to  the  hysteresis  at  T  =  263°  K  but  with  less  accuracy  at  the  higher  temperatures. 


Example  2  (PMN): 

To  provide  a  second  illustration  of  the  model  attributes  over  a  more  extensive  temperature  range,  we 
consider  its  capability  for  predicting  the  dielectric  behavior  of  PMN  data  collected  by  Glazounov  et  al.  [7] . 
The  data  set  under  consideration  is  plotted  in  Figure  9.  For  temperatures  ranging  from  T  =  133°  K  to 
T  =  313°  K,  the  input  fields  were  prescribed  to  have  peak  values  of  3.52  MV/rn  and  frequencies  of  0.1  Hz 
to  minimize  frequency  effects.  The  data  exhibits  significant  hysteresis  in  the  low  temperature  ferroelectric 
regimes  with  decreasing  hysteresis  and  saturation  values  as  the  temperatures  increase. 

To  illustrate  the  thermal  dependence  of  the  hysteresis  loss  parameter  k(T)  for  this  temperature  range, 
the  coercive  fields  are  plotted  as  a  function  of  temperature  in  Figure  10.  Additionally,  the  best  least  squares 
linear  and  quadratic  fits  to  the  data  are  plotted  to  motivate  the  choice  of  p  in  the  relation 

k{T)  =  k(Tf  —  T)v  , 

where  k(T)  «  EC(T),  derived  in  Section  3.2.  It  is  observed  that  for  temperatures  above  213°  K,  the  coercive 
data  exhibits  a  nearly  linear  dependence  on  temperature  and  one  can  employ  p  =  1  as  in  Example  1.  For 
the  full  temperature  range,  however,  a  linear  model  will  significantly  underestimate  the  coercive  field  at  low 
temperatures  and  a  quadratic  model  is  required  for  accurate  predictions;  hence  for  the  model  development, 
we  employed  p  =  2.  Finally,  the  data  in  Figure  10  indicates  that  the  critical  freezing  temperature  for  this 
material  lies  between  233°  K  and  253°  K,  and  we  employed  Tf  =  242°  K  in  the  model. 

The  bulk  Curie  temperature  was  taken  to  be  Tcur  =  398°  K.  The  remaining  model  parameters  were 
estimated  through  a  least  squares^  fit  to  the  data.  The  model  predictions  obtained  with  the  parameter 
set  a  —  1.8  x  108,a  =  2.5  x  108,fc  =  190,  Ps  =  0.38, /?  =  1.42  x  10-8  and  c  =  0.1  are  compared  to  the 
PMN  data  in  Figure  9.  It  is  observed  that  the  model  quantifies  both  the  hysteretic  dielectric  behavior 
in  the  ferroelectric  regime  and  the  decaying  anhysteretic  response  above  the  freezing  temperature.  The 
primary  discrepancy  between  the  data  and  model  predictions  lies  in  the  elbow  of  the  hysteresis  curve 
where  the  model  predicts  a  more  gradual  transition  from  saturation  to  coercivity  than  is  observed  in  the 
data.  It  is  hypothesized  that  this  is  due  to  domain  rotation  mechanisms  which  are  not  incorporated  in  the 
anhysteretic  model.  We  note  that  the  model  fits  at  any  given  temperature  can  be  improved  by  optimizing 
the  parameters  for  that  temperature.  However,  the  specified  parameter  set  provides  a  good  fit  throughout 
the  temperature  range,  with  thermal  dependence  in  a(T),  a(T),  k(T )  and  PS(T)  incorporated  solely  through 
the  expressions  (37). 
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Figure  10.  Coercive  PMN  data  from  Glazounov  et  al.  [7]  plotted  as  a  function  of  temperature  and  the 
best  least  squares  linear  and  quadratic  fits  to  the  data. 

5  Concluding  Remarks 

This  paper  addresses  the  modeling  of  temperature-dependent  constitutive  nonlinearities  and  hysteresis  in 
relaxor  ferroelectric  materials.  As  in  the  fixed-temperature  case,  the  model  is  constructed  in  two  steps: 
(i)  the  modeling  of  the  anhysteretic  dielectric  behavior  through  Boltzmann  and  thermodynamic  principles, 
and  (ii)  the  modeling  of  hysteresis  through  the  quantification  of  energy  required  to  bend  and  translated 
domain  walls  pinned  at  inclusions  inherent  to  the  materials.  The  microscopic  anhysteretic  model  is  based 
on  the  assumption  that  the  materials  are  comprised  of  micropolar  regions  having  a  broad  spectrum  of 
Curie  temperature  due  to  their  chemical  heterogeneity.  Consideration  of  the  free  energy  in  the  presence 
of  variations  in  cation  distributions  yields  an  Ising  spin  model  representing  the  thermodynamics  of  the 
microregions.  Based  on  a  “random  layer”  model  for  B-site  orientations,  this  then  provided  a  mechanism  for 
computing  the  saturation  polarization  and  distribution  of  micropolar  regions  as  a  function  of  temperature. 
A  nonlinear  model  characterizing  the  aggregate  material  behavior,  in  the  absence  of  hysteresis,  was  then 
obtained  by  combining  the  micropolar  model  with  macroscopic  constitutive  relations.  To  quantify  the 
hysteresis  inherent  to  the  materials  at  low  temperatures,  the  energy  required  to  reorient  dipoles,  and 
hence  move  domain  walls,  was  quantified  as  a  function  of  temperature.  The  resulting  combined  model 
characterizes  the  hysteresis  and  constitutive  nonlinearities  exhibited  by  the  materials  from  low  temperature 
ferroelectric  regimes  through  the  superparaelectric  phase  where  the  materials  exhibit  saturation  effects  but 
negligible  hysteresis. 

The  predictive  capabilities  of  the  model  were  illustrated  through  comparison  with  PMN  and  PMN-PT- 
BT  data  collected  under  quasistatic  conditions.  For  both  data  sets,  the  model  quantifies  the  hysteresis 
observed  at  low  temperatures  as  well  as  the  decrease  in  hysteresis  which  occurs  as  temperatures  increase. 
It  also  quantifies  the  transition  to  anhysteretic  behavior  and  the  reduction  in  saturation  polarization  which 
occurs  as  temperatures  increase.  The  full  temperature  and  field  dependence  is  incorporated  in  the  model 
and  the  differing  model  responses  were  obtained  with  one  set  of  parameters.  The  sole  inputs  to  the  model 
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are  the  input  field  and  the  operating  temperature.  The  small  number  of  parameters  (7)  and  their  asymptotic 
relation  to  measurable  quantifies  (e.g.,  coercive  field  Ec.  slope  at  field  reversal)  facilitates  the  construction 
of  the  model  as  well  as  the  updating  of  parameters  to  accommodate  changing  operating  conditions. 

Due  to  the  focus  of  the  paper  on  the  quantification  of  the  thermal  dependence  on  the  dielectric  prop¬ 
erties  of  relaxor  materials,  the  examples  were  chosen  to  have  fixed  maximum  field  inputs  and  variable 
temperatures.  The  capability  of  the  model  to  accommodate  variable  field  inputs  is  documented  for  various 
PZT  compounds  in  [21]  through  a  comprehensive  set  of  examples. 

We  note  that  the  model  developed  here  is  quasistatic  in  nature  and  hence  does  not  yet  incorporate 
the  physical  mechanisms  which  produce  certain  frequency  effects  in  relaxor  materials.  Care  should  thus 
be  exhibited  if  the  model  is  employed  in  variable  or  high  frequency  regimes.  While  certain  choices  for 
parameters  may  produce  accurate  model  fits  at  fixed  high  frequencies,  this  should  be  attributed  to  the 
mathematical  flexibility  of  the  model  rather  than  incorporated  physics.  The  development  of  a  commensu¬ 
rate  physics-based  dynamical  model  is  under  current  investigation. 

Finally,  the  ODE  nature  of  the  model  makes  it  amenable  to  inversion  through  the  consideration  of  a 
complementary  ODE  in  a  manner  analogous  to  that  described  in  [17].  This  facilitates  the  construction 
of  an  inverse  compensator  which  can  be  used  for  linear  control  design  [24,  27].  The  application  of  these 
techniques  for  linear  control  implementation  for  piezoceramic  actuators  which  exhibit  hysteresis  is  also 
under  investigation. 
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